Accuracy of density functionals for molecular electronics: the Anderson junction 
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The exact ground-state exchange-correlation functional of Kohn-Sham density functional theory 
yields the exact transmission through an Anderson junction at zero bias and temperature. The 
exact impurity charge susceptibility is used to construct the exact exchange-correlation potential. 
We analyze the successes and limitations of various types of approximations, including smooth and 
discontinuous functionals of the occupation, as well as symmetry-broken approaches. 



Since the pioneering experiments of Reed and Tour 
on dithiolated benzene [TJ, there has been tremendous 
progress in the ability to both create and characterize 
[2] organic molecular junctions. But accurate simulation 
of these devices remains a challenge, both theoretically 
and computationally !3|. The essential physics has been 
well understood since the ground-breaking work of Lan- 
dauer and Biittiker [4] [5] in the context of mesoscopic 
devices, including both Coulomb blockade and Kondo 
effects [SI [7J . Calculations with simple model Hamilto- 
nians demonstrate such effects at a qualitative level [5]. 
On the other hand, organic molecules connected to metal 
leads [9] require hundreds of atoms and thousands of ba- 
sis functions for a sufficiently accurate calculation of their 
total energy, geometry, and single-particle states. Such 
conditions are routine for modern density functional the- 
ory (DFT) calculations [10 , but the ability of present 
functional approximations to predict accurate currents 
remains an open question |11) . 

The standard DFT method for calculating current 
through such a device is to perform a ground-state Kohn- 
Sham (KS) DFT calculation [12] on a system upon which 
a difference between the chemical potentials of the left 
and right leads has been imposed (the applied bias), and 
calculate the transmission through the KS potential us- 
ing the Landauer-Biittiker formula. But there is nothing 
in the basic theorems of DFT that directly implies that 
such a calculation would yield the correct current, even 
if the exact ground-state functional were used. 

The limit of weak bias is more easily analyzed than the 
general case, because the Kubo linear response formalism 
applies [T31 HI] . In that case one finds that, in princi- 
ple, there are exchange-correlation (XC) corrections to 
the current in the standard approach [15) . but little is 
known about their magnitude |16[ 117) . Even without 
these corrections, one can ask if the standard approx- 
imations used in most ground-state DFT calculations 
(i.e., generalized gradient approximations [TJ] and hy- 
brids of these with Hartree-Fock exchange [HHHD]) are 
sufficiently accurate for transport purposes. The answer 
appears definitively no! Because of self-interaction errors, 
such approximations are well-known [21] to produce po- 



tentials with incorrectly positioned KS eigenvalues, both 
occupied and unoccupied. These errors become severe 
when the molecule is only weakly coupled to the leads 
[2"2] . Calculated transmission can be too large by several 
orders of magnitude due to this incorrect positioning of 
the levels. Recent calculations [23] using beyond-DFT 
techniques to correctly position the levels show greatly 
improved agreement with experiment. 

But this progress returns us to the earlier concern: 
Even with an exact ground-state XC functional, are there 
XC corrections to the Landauer-Biittiker result? The an- 
swer appears to be yes in general [15] . but in a previous 
work |24j we argued that, under a broad range of condi- 
tions applicable to typical experiments, such corrections 
can vanish. This result was shown by exact calculations 
on an impurity model (Anderson model) employing the 
exact XC functional. In the present work, we analyze dif- 
ferent approximate treatments, applied to the Anderson 
junction, and calculate their errors. The implications for 
DFT calculations of transport in general are discussed. 

The Anderson model %ZE\ is a single interacting site 
(C) connected to two non-interacting electrodes (L,R). 
The Hamiltonian of the system is T-L = He + Ht + %l,r- 
Each lead is represented by a non-interacting Fermi gas: 
"Hl,r = J2kaeL r E ko-n a , with chemical potential \i and 
the central interacting site is: He — £ (?V + n±) + Unfn±, 
where fi a = d\d a is the number operator for spin a and 
U is the charging energy representing on-site interaction. 
Ht is the tunneling between leads and the central site. 
The tunneling width T is a constant in the broad-band 
limit. A schematic is shown in Fig. [T] Real molecules 
can be mapped onto the Anderson model [2H1 [57] . 

In a previous work [53] , we calculated the exact relation 
between occupancy on the central site and on-site energy 
e for an Anderson junction, using the Bethe ansatz (BA) 
[25] . We showed that exact KS DFT yields the exact 
transport at zero temperature and in the linear response 
regime, although the KS spectral function differs from 
the exact one away from the Fermi energy. This is be- 
cause the Anderson junction has only one site and trans- 
mission is a function of occupation number due to the 
Friedel-Langreth sum rule [2^1 ED] ■ Thus, for this sim- 
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FIG. 1. A cartoon for Anderson model. The model consists 
of two featureless leads and a central region with on-site in- 
teraction U. F is the tunneling width. Two many-body levels 
of the central region are shown. 
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pie model, all failures of approximate XC calculations of 
transmission can be attributed to failures to reproduce 
the exact occupation number, i.e., there are no XC cor- 
rections to the standard practice of applying KS DFT to 
the ground-state and finding transmission through the 
single-particle potential. On the other hand, the stan- 
dard approximations in use in DFT calculations of trans- 
port have a variety of shortcomings. The most prominent 
one, as we shall see, is the lack of a discontinuity in the 
XC potential with particle number [3T] . 

Before studying approximations, we refine our previous 
numerical fit of BA results, using analytic results from 
many-body theory. We re-introduce [35] reduced vari- 
ables y = T /U , which measures the ratio of lead-coupling 
to the onsite Coulomb repulsion, while x = (fi — e)/U is 
the difference between the leads' chemical potential and 
the onsite level energy, in units of U. For x < 0, the cen- 
tral site is above the chemical potential, at x = they 
match. 

The occupation in the KS system is given by self- 
consistent solution of the KS equation for occupation: 
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where the KS level is written as 
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with the second term being the Hartree contribution and 
the third being the XC contribution (in fact, only cor- 
relation, as exchange is zero for this model), which is a 
function of the occupation. Considered in reverse, this is 
a definition of the exact e X c > if the occupation is known, 
as it is from the BA solution. The KS transmission is 
then 
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and matches the true transmission in the many-body 
system, by virtue of the sum-rule. The exact ground- 
state functional yields the exact transmission, including 



FIG. 2. Dimensionless susceptibility x c = Uxc(l) as a func- 
tion of y = U/T for the Anderson junction [exact, [5,6]-Pade 
fit (see text), RHF and UHF]. 



the Kondo plateau at zero temperature and weak bias 

EH 132] . 

As shown in Ref. [21], the XC potential can be very 
accurately parametrized with the form: 
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The tan -1 term jumps by 7r as (n c ) passes through 1, 
leading to discontinuous behavior with occupation. Thus 
a determines the width of this region, while a determines 
its strength. Both a and a are functions of y = T/U and 
were extracted numerically by fitting to the exact solu- 
tion, and were roughly fit by simple Pade approximations 
there. 

However, we can greatly improve the fit of a. A central 
object in the Anderson junction is the charge susceptibil- 
ity, Xc((^c)) — d (n c ) /d/i. At half-filling, this is known 
analytically [331 [34]: 
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where \ c = UxcO-) is dimensionless, and is plotted in 
Fig. [2] This curve can be readily fit to a [5,6] Pade form: 



xr d (y) = E a ^/E 6 ^ fc ' 
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whose 11 independent coefficients are chosen to recover 
the Taylor-expansion around y — (strongly-correlated 
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limit exactly to 5 orders, around y — > oo to 6 or- 
ders, and are given in Table |TJ The weak-correlation limit 
can also be extracted via the Yosida-Yamada perturba- 
tive approach |37H39| . The quantity \ c has the physical 
meaning of the slope at particle-hole symmetry point in 



the (n c ) vs. (fi — e)/U curve (see Figs. [3] and [4]). It has 
a maximum at about y = 0.291, and as y varies from oo 
(weakly-correlated limit) to (strongly-correlated limit), 
the slope at the symmetric point increases at first. Be- 
yond the maximum value, the slope decreases, and the 
Coulomb blockade pleateau gradually develops. 



TABLE I. Coefficients in the [5,6]-Pade approximation [Eq. jfjj]. 
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By taking derivatives on both sides of Eq. Q , the two 
coefficients a and a are constrained by \ c : 

g = . TV- (7) 

(2/xc - yn + a — l) 

Retaining the simple form of Ref. [21], a [0,1] Pade, a = 
1/(1 + 5.68y), we determine a from the Pade fit to the 
susceptibility and Eq. Q. This yields highly accurate 
occupations, KS potentials, and transmissions, including 
the Kondo plateau. It agrees very well with the numerical 
fit to the BA results of Ref. [24], and matches more 
closely than the simpler analytic fit used there. 

We now move on to the central topic of this work, 
which is the accuracy of approximate functional treat- 
ments. In such treatments, e X c is approximated as a 
function of (n c ) in Eq. (J2J, and the resulting Eq. |l]) is 
solved self-consistently for (n c ). The simplest such ap- 
proximation is to simply set e xc = 0, i.e., Hartree-Fock 
(HF), and should be accurate when correlation is weak. 
In Fig. [3] we plot several quantities for U = V, both ex- 
actly and in HF, showing that HF is very accurate here. 
We find [25]: 
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1 In Refs. I35| and I36| . this expansion was reported incorrectly, 
with minus sign on the second term. We believe it is corre- 
sponding to Wilson ratio R = 1, however this is not true in 
strongly-correlated limit. 
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FIG. 3. Upper panel: transmission as a function of x = 
(fj, — e)/U ; middle panel: occupation as a function of x; lower 
panel: KS potential as a function of occupation. Results are 
shown for Bethe ansatz or exact KS DFT (exact), Hatree- 
Fock (HF), and discontinuous approximation [disc, Eq. ( 10 1] . 
U — r in all cases. 



where 7 = 3 — 7r 2 /4 exactly, but 7 = 1 in HF. Thus 
we regard U < T as the weakly correlated regime. On 
the other hand, in Fig. [4] we show the same plots for 
U = 10 r. Now, in the exact occupation, the slope near 
x = 0.5 is much weaker, leading to a transmission plateau 



4 




FIG. 4. Upper panel: transmission as a function of x = 
(/Lt— t)/U; middle panel: occupation as a function of a;; lower 
panel: KS potential as a function of occupation. Results are 
shown for Bethe ansatz or exact KS DFT (exact), restricted 
Hatree-Fock (RHF), unrestricted Hartree-Fock (UHF), and 
discontinuous approximation [disc, Eq. (|10[)]. U = 10r in all 
cases. 



(the Kondo plateau) for < x < 1. The plateau effect 
is missed entirely by HF, because of the too-smooth de- 
pendence (in fact, linear) of its KS level on occupation 
(see bottom panel). Note that at temperatures equal 
to or above the Kondo temperature, the Kondo effect 
is destroyed, and the central plateau in transmission is 
replaced by two Hubbard peaks around x = and 1. 
Then the behavior of the HF curve is exactly as qualita- 
tively predicted in Ref. [15] . smearing out the two sharp 
features into one peak midway between them. This is be- 
cause the KS level shifts linearly with occupation in HF, 
instead of more suddenly with occupation in the exact 
solution. More generally, all smooth density function- 
als, such as the local density approximation |12| and the 
generalized gradient approximation |18) . suffer from the 
same qualitative failure, and so would produce incorrect 
peaks centered at x = 0.5. All these errors arise from the 
approximations to the functional; the exact ground-state 
functional reproduces the exact occupation by construc- 
tion, and so yields the exact transmission. 

There have thus been several suggestions [40] to incor- 
porate the discontinuous behavior with occupation into 
approximations in transport calculations. At the prac- 
tical level, Toher et al. [22] showed in a model calcula- 
tion how self-interaction corrections would greatly sup- 
press zero-bias conductance in local density approxima- 
tion calculations for molecules weakly coupled to leads. 
More recently, Bethe Ansatz Local Density Approxima- 



tion (BALDA) 40J and variations [H] have been used 
to impose discontinuous behavior on the levels. For sim- 
ple models, all of these can be considered as LDA+L 7 - 
like. The methodology of LDA+C/ [43] has become in- 
creasingly popular in recent years, especially for those 
focused on moderately correlated systems such as tran- 
sition metal oxides, for which LDA and GGA often have 
zero KS band gap. In some fashion, a Hubbard U is 
added to some orbitals of a DFT Hamiltonian. Some- 
times U is regarded as an empirical parameter, while 
others have found self-consistent prescriptions. In any 
event, despite not fitting in the strict DFT framework, it 
is a method borne of practical necessity for many situa- 
tions [43] . 

To gain a qualitative understanding of the effects of 
such models, we define a very simple XC potential that 
has a discontinuity. To do this, we simply take the 
Hartree form, symmetrize it around the half- filled point, 
and replace U by a screened U. We find that a sim- 
ple fit U — U/(l + 0.25/y) works well. U being different 
from U and particle-hole symmetry guarantee an explicit 
derivative discontinuity of e s with respect to occupation 
number. This yields 



e s [n] = ~Un6{l-n) + 



U+ l -U{n-2) 



9(n-l), (10) 



where 0(x) is the Heaviside theta function and for sim- 
plicity, n is just (n c ). 

While this model does contain a discontinuity, and 
yields the exact result as y — > 0, curing the worst de- 
fects of HF, it misses entirely the finite slope of the KS 
potential at half-filling for finite U, which is determined 
by the susceptibility. The explicit derivative discontinu- 
ity is exact in the strongly-correlated limit with infinite 
U/T, but should be "rounded" in finite U/T [Mid], or 
in finite temperature [44] . To see this for finite (but very 
large) U/T, in Fig. [5) we show similar results as in Figs. 
[3]and|4j but with U = 100r and we only show the region 
around (n c ) = 1 at x = 0, where the rounded derivative 
discontinuity occurs. The transmission is accurate both 
for weak and strong correlation, but is not so everywhere 
in between. In particular, it is overestimated for (n c ) just 
above (and just below 1) for U = 10r because of this 
lack of a finite slope. This is where we expect the great- 
est errors in such models, but the region of inaccuracy 
(on the scale of x) shrinks as U/T —> oo. 

Finally, we discuss a different class of approximations. 
A well-known (and much debated) technique for mim- 
icking strong correlation is to allow a mean-field calcula- 
tion to break symmetries that are preserved in the exact 
calculation. Perhaps the most celebrated prototype of 
such a calculation is for HF applied to an H2 molecule 
with a large bond distance. At a crucial value of the 
bond distance (called the Coulson-Fischer point), an un- 
restricted calculation, i.e., one that allows a difference 
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in spin occupations, yields a lower energy than the re- 
stricted one. This remains the case for all larger separa- 
tions, and the unrestricted solution correctly yields the 
sum of atomic energies as R — > oo, whereas the restricted 
Hartree-Fock (RHF) solution dissociates to unpolarized 
H atoms with the wrong energies. This is the celebrated 
symmetry dilemma: with a mean-field approximation, 
for large separations, one can either get the right symme- 
try (RHF) or the right energy [unrestricted Hartree-Fock 
(UHF)], but not both. The same issues arise in approx- 
imate DFT treatments of this problem [33] ■ Of course, 
the exact functional manages to get the correct energy 
with the correct symmetry, and there have been many 
attempts to reproduce this with various more sophisti- 
cated approximations. But a more pragmatic approach 
is to accept the results as they are, interpreting the good 
energetics as the result of applying the approximate func- 
tional to a frozen fluctuation of the system. The true 
ground-state wavefunction fluctuates between configura- 
tions with one spin and then the other (left and right 
localized for stretched H2), and the true ground-state 
density has unbroken symmetry. But the approximate 
functionals give most accurate energies when applied to 
the frozen fluctuations. Thus, one can interpret both the 
total density and energy as being accurate from such a 
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FIG. 5. Upper panel: transmission as a function of x = 
(fj.— e)/U; middle panel: occupation as a function of x; lower 
panel: KS potential as a function of occupation. Results are 
shown for Bethe ansatz or exact KS DFT (exact), unrestricted 
Hartree-Fock (UHF), and discontinuous approximation [disc, 
Eq. ( 10 1] . Also shown in the upper panel is UHF results for 
transmission using (incorrect) spin densities [UHF(SB), with 
symmetry breaking], and (m) = (rif) — (nj.) for UHF as a 
function of x as an inset in the middle panel. U — 100F in all 
cases, and only the region near (no) = 1 and x = is shown. 



calculation, but not the individual spin-densities. In fact, 
an alternative approach is to interpret another variable, 
such as the ontop pair density, as being accurately ap- 
proximated in such treatments |45j. 

We apply the same reasoning to the Anderson junction, 
just as was done by Anderson when creating the model 
we are using |25| . The symmetries are different, but the 
principle is the same. We allow the mean-field calculation 
to break spin-symmetry if this leads to lower energy on 
the central site, with spin equations: 



(nt) 
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■ arctan £ - U ^ ' ) 

(11) 

and reverse for (fij,), and (n c ) = (n-j-) + (n±). Again, the 
simplest calculation is UHF, where e xc = 0. The solu- 
tions are identical to those found in the original problem 
by Anderson [35]. For y > l/V, i.e., U < irT, there is no 
spontaneous symmetry-breaking, and UHF=RHF. But 
beyond that critical value, the spin-density difference be- 
comes finite, and the unrestricted solution differs. Define 
the density difference as (to) = (n-j-) — (n±) in UHF, which 
satisfies: 



tan — (to) = 



2y 



UHF 



(12) 



For y > l/V, (to) = 0, but otherwise a solution with (m) 
finite exists. In all cases, we take only the total density 
from the UHF calculation, and we know the true (to) = 
always. In particular, as y — ¥ (strong correlation), 
Xc — > with the correct linear term: 



Xc^ (8/7r)y + (96 7 /tt V 



(13) 



where 7 = 1 in the exact solution, but 7 = 1/3 in UHF. 
So UHF recovers the leading term. The green curve in 
Fig. [2]shows the UHF value of % c , demonstrating both its 
accuracy for both strong and weakly correlated systems, 
and the discontinuous change at 1/7T. 

Even beyond the "Coulson-Fisher point" of 1/tt, the 
symmetry-breaking only occurs for < (to) < 1, i.e., 
outside this region, the UHF solution is that of RHF, as 
can be seen in the inset of middle panel in Fig. [5j But the 
density is very accurately given by UHF (considering the 
scale of horizontal axis), and the KS potential develops 
the correct derivative discontinuity as y — > 0. 

To demonstrate this accuracy, we plot the correspond- 
ing transmissions in Fig. [4j using Eq. ^. The fig- 
ure shows how the transmission using (n c ) from UHF is 
almost exact (considering the scale of horizontal axis). 
To demonstrate the error in ignoring the fact that the 
UHF produces incorrect spin densities, we also plot the 
transmission through such a solution, which is completely 
wrong (see dark red curve in upper panel of Fig. [5j only 
one peak is present because only region near x = and 
(n c ) = 1 is shown there). Our results are consistent with 
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those of [25], justifying the use of the broken symmetry 
solution to deal with strong correlation. 

To summarize, we have studied approximate treat- 
ments of the zero-temperature weak-bias conductance 
of the Anderson junction. RHF and approximate DFT 
treatments work well for weak correlation, but fail for 
moderate and strong correlation because of the smooth 
dependence of their KS potentials on occupation num- 
bers. Imposing an explicit discontinuity consistent with 
particle-hole symmetry can yield a discontinuity with oc- 
cupation which guarantees correct behavior in the strong 
correlation limit. This also greatly improves results for 
moderate correlation, but still contains errors. Finally, 
simple symmetry-breaking in UHF produces remarkably 
accurate conductances, once the transmission is calcu- 
lated as if the symmetry had not been broken. 
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